# Импротируем необходимые функции
import numpy as np # работа с массивами и линейной алгеброй
import matplotlib.pyplot as plt # для отрисовки графиков
import pandas as pd # для чтения и работы с данными
from scipy.optimize import curve_fit # фитирующая процедура

def psi_f(x):
    return np.exp(-x**2/2)/np.sqrt(np.sqrt(np.pi))

def integrate_hist(hist, bin_edges):
    integral = 0
    for i in range(len(hist)):
        integral += hist[i]*(bin_edges[i+1] - bin_edges[i])
    return integral

def distr(df):
    pedestrians, bin_edges = np.histogram(df, bins = 1000)
    #print bin_edges
    
    width = bin_edges[1] - bin_edges[0]
    plt.figure(figsize=(7,7))
    plt.bar(bin_edges[:-1], pedestrians, width = width)
    plt.xlabel('Ось x'.decode('utf-8'))
    plt.ylabel('Количество пешеходов (N)'.decode('utf-8'))
    plt.title('Распределение пешеходов по оси'.decode('utf-8'))
    plt.legend()
    plt.show()
    
    pedestrians_sq = pedestrians*pedestrians
    norm = integrate_hist(pedestrians_sq, bin_edges)
    psi_sq = pedestrians_sq/norm
    psi = pedestrians/np.sqrt(norm)

    plt.figure(figsize=(7,7))
    plt.bar(bin_edges[:-1], psi, width = width)
    plt.xlabel('x'.decode('utf-8'))
    plt.ylabel('$\psi(x)$'.decode('utf-8'))
    plt.title(''.decode('utf-8'))
    ###
    x = np.arange(-4, 4, 0.1)
    plt.plot(x, psi_f(x), color = 'black', label = '$\psi(x)$, теоретическая'.decode('utf-8'))
    plt.legend()
    plt.show()
    
    plt.figure(figsize=(7,7))
    plt.bar(bin_edges[:-1], psi_sq, width = width)
    plt.xlabel('x'.decode('utf-8'))
    plt.ylabel('$\psi^2(x)$'.decode('utf-8'))
    plt.title('Плотность вероятности'.decode('utf-8'))
    plt.plot(x, psi_f(x)**2, color = 'black', label = '$\psi^2(x)$, теоретическая'.decode('utf-8'))
    plt.legend()
    plt.show()
    
def drawNE(df):
    plt.figure(figsize=(20,7))
    plt.scatter(np.arange(df.shape[0]), df['N'], linewidths=0.1, label = 'Количество пешеходов'.decode('utf-8') )
    plt.xlabel('Номер итерации'.decode('utf-8'))
    plt.ylabel('Количество пешеходов (N)'.decode('utf-8'))
    plt.title('Зависимость N от итерации'.decode('utf-8'))
    plt.legend()
    plt.grid()
    plt.show()
    
    plt.figure(figsize=(20,8))
    plt.plot(np.arange(df.shape[0]), df['E'], label = 'Данные'.decode('utf-8'))
    plt.xlabel('Номер итерации'.decode('utf-8'))
    plt.ylabel('E'.decode('utf-8'))
    plt.title('Зависимость энергии от итерации'.decode('utf-8'))
    plt.legend()
    plt.grid()
    plt.show()
    
    plt.figure(figsize=(20,8))
    plt.plot(np.arange(df.shape[0]), df['E_r'], label = 'Данные'.decode('utf-8'))
    plt.xlabel('Номер итерации'.decode('utf-8'))
    plt.ylabel('E_r'.decode('utf-8'))
    plt.title('Зависимость опорного потенциала от итерации'.decode('utf-8'))
    plt.legend()
    plt.grid()
    plt.show()
    
E_ds = pd.DataFrame(columns=['ds', 'E'])

def analyse(frame, frameNE, ds, E_ds, begin = 500):
    distr(frame)
    drawNE(frameNE)
    E = frameNE[begin:]['E'].mean()
    s_E = np.std(frameNE[begin:]['E'])/np.sqrt(frameNE[begin:].shape[0])
    E_r = frameNE[begin:]['E_r'].mean()
    s_E_r = np.std(frameNE[begin:]['E_r'])/np.sqrt(frameNE[begin:].shape[0])
    print E, s_E
    print E_r, s_E_r
    s = pd.Series([ds,E],index=['ds','E'])
    E_ds = E_ds.append(s,ignore_index=True)
    return E_ds

frame = pd.read_csv('D:/Backup(work)/data_output/output(1000000x1000)_ds=0.11/output.csv')
frameNE = pd.read_csv('D:/Backup(work)/data_output/output(1000000x1000)_ds=0.11/outputNE.csv')
E_ds = analyse(frame, frameNE, ds = 0.11, E_ds = E_ds)
del frame
del frameNE

0.495383074 4.61381833306e-05
0.500871786 0.000277473401522

frame = pd.read_csv('D:/Backup(work)/data_output/output(1000000x1000)_ds=0.12/output.csv')
frameNE = pd.read_csv('D:/Backup(work)/data_output/output(1000000x1000)_ds=0.12/outputNE.csv')
E_ds = analyse(frame, frameNE, ds = 0.12, E_ds = E_ds)
del frame
del frameNE

0.495062006 2.30052662651e-05
0.501981934 0.000251843564816

E_ds

	ds	E
0	0.11	0.495383
1	0.12	0.495062

frame = pd.read_csv('D:/Backup(work)/data_output/output(1000000x1000)_ds=0.15/output.csv')
frameNE = pd.read_csv('D:/Backup(work)/data_output/output(1000000x1000)_ds=0.15/outputNE.csv')
E_ds = analyse(frame, frameNE, ds = 0.15, E_ds = E_ds)
del frame
del frameNE

0.492044908 3.57702241406e-05
0.502984848 0.000191987294918

frame = pd.read_csv('D:/Backup(work)/data_output/output(1000000x1000)_ds=0.095/output.csv')
frameNE = pd.read_csv('D:/Backup(work)/data_output/output(1000000x1000)_ds=0.095/outputNE.csv')
E_ds = analyse(frame, frameNE, ds = 0.095, E_ds = E_ds)
del frame
del frameNE

0.495843906 3.17241787022e-05
0.499505012 0.000310959743928

frame = pd.read_csv('D:/Backup(work)/data_output/output(1000000x1000)_ds=0.096/output.csv')
frameNE = pd.read_csv('D:/Backup(work)/data_output/output(1000000x1000)_ds=0.096/outputNE.csv')
E_ds = analyse(frame, frameNE, ds = 0.096, E_ds = E_ds)
del frame
del frameNE

0.49616667 2.95552506706e-05
0.50017879 0.000311246394401

frame = pd.read_csv('D:/Backup(work)/data_output/output(1000000x1000)_ds=0.097/output.csv')
frameNE = pd.read_csv('D:/Backup(work)/data_output/output(1000000x1000)_ds=0.097/outputNE.csv')
E_ds = analyse(frame, frameNE, ds = 0.097, E_ds = E_ds)
del frame
del frameNE

0.496018188 2.73999552064e-05
0.499670628 0.000308254392571

frame = pd.read_csv('D:/Backup(work)/data_output/output(1000000x1000)_ds=0.098/output.csv')
frameNE = pd.read_csv('D:/Backup(work)/data_output/output(1000000x1000)_ds=0.098/outputNE.csv')
E_ds = analyse(frame, frameNE, ds = 0.098, E_ds = E_ds)
del frame
del frameNE

0.496273366 3.51456105949e-05
0.500470488 0.000321138522092

frame = pd.read_csv('D:/Backup(work)/data_output/output(1000000x1000)_ds=0.099/output.csv')
frameNE = pd.read_csv('D:/Backup(work)/data_output/output(1000000x1000)_ds=0.099/outputNE.csv')
E_ds = analyse(frame, frameNE, ds = 0.099, E_ds = E_ds)
del frame
del frameNE

0.496013276 2.27433007202e-05
0.499874848 0.000311800388739

frame = pd.read_csv('D:/Backup(work)/data_output/output(1000000x1000)_ds=0.08/output.csv')
frameNE = pd.read_csv('D:/Backup(work)/data_output/output(1000000x1000)_ds=0.08/outputNE.csv')
E_ds = analyse(frame, frameNE, ds = 0.08, E_ds = E_ds, begin = 800)
del frame
del frameNE

0.497641675 2.61326777594e-05
0.499825855 0.000625934763909

frame = pd.read_csv('D:/Backup(work)/data_output/output(1000000x1000)_ds=0.07/output.csv')
frameNE = pd.read_csv('D:/Backup(work)/data_output/output(1000000x1000)_ds=0.07/outputNE.csv')
E_ds = analyse(frame, frameNE, ds = 0.07, E_ds = E_ds, begin = 1000)
del frame
del frameNE

0.497182122 2.05621986721e-05
0.49849758 0.000461492505594

frame = pd.read_csv('D:/Backup(work)/data_output/output(1000000x1000)_ds=0.06/output.csv')
frameNE = pd.read_csv('D:/Backup(work)/data_output/output(1000000x1000)_ds=0.06/outputNE.csv')
E_ds = analyse(frame, frameNE, ds = 0.06, E_ds = E_ds, begin = 1500)
del frame
del frameNE

0.49723906 3.32053792148e-05
0.49722688 0.000508842110988

frame = pd.read_csv('D:/Backup(work)/data_output/output(1000000x2000)_ds=0.05/output.csv')
frameNE = pd.read_csv('D:/Backup(work)/data_output/output(1000000x2000)_ds=0.05/outputNE.csv')
E_ds = analyse(frame, frameNE, ds = 0.05, E_ds = E_ds, begin = 1500)
del frame
del frameNE

0.496872982 2.01907748081e-05
0.496243498 0.000618111813623

frame = pd.read_csv('D:/Backup(work)/data_output/output(1000000x4000)_ds=0.04/output.csv')
frameNE = pd.read_csv('D:/Backup(work)/data_output/output(1000000x4000)_ds=0.04/outputNE.csv')

E_ds = analyse(frame, frameNE, ds = 0.04, E_ds = E_ds, begin = 3000)
del frame
del frameNE

0.496631159 2.36421893174e-05
0.493893758 0.000552316819061

frame = pd.read_csv('D:/Backup(work)/data_output/output(1000000x4000)_ds=0.03/output.csv')
frameNE = pd.read_csv('D:/Backup(work)/data_output/output(1000000x4000)_ds=0.03/outputNE.csv')
E_ds = analyse(frame, frameNE, ds = 0.03, E_ds = E_ds, begin = 3500)
del frame
del frameNE

0.493640188 9.9354569755e-06
0.489054732 0.000951642280015

frame = pd.read_csv('D:/Backup(work)/data_output/output(1000000x1500)_ds=0.065/output.csv')
frameNE = pd.read_csv('D:/Backup(work)/data_output/output(1000000x1500)_ds=0.065/outputNE.csv')
E_ds = analyse(frame, frameNE, ds = 0.065, E_ds = E_ds, begin = 1000)
del frame
del frameNE

0.497278778 1.84619529149e-05
0.497968004 0.000490557220165

frame = pd.read_csv('D:/Backup(work)/data_output/output(1000000x1500)_ds=0.075/output.csv')
frameNE = pd.read_csv('D:/Backup(work)/data_output/output(1000000x1500)_ds=0.075/outputNE.csv')
E_ds = analyse(frame, frameNE, ds = 0.075, E_ds = E_ds, begin = 1000)
del frame
del frameNE

0.496998904 2.49333410029e-05
0.4986877 0.000411465900405

frame = pd.read_csv('D:/Backup(work)/data_output/output(1000000x1500)_ds=0.085/output.csv')
frameNE = pd.read_csv('D:/Backup(work)/data_output/output(1000000x1500)_ds=0.085/outputNE.csv')
E_ds = analyse(frame, frameNE, ds = 0.085, E_ds = E_ds, begin = 1000)
del frame
del frameNE

0.496556478 1.67013552454e-05
0.498796588 0.000366256550757

E_ds = E_ds.sort_values(by=['ds'])
E_ds.to_csv('E_ds.csv')
plt.plot(E_ds['ds'], E_ds['E'])
plt.grid()
plt.show()

E_ds

	ds	E
13	0.030	0.493640
12	0.040	0.496631
11	0.050	0.496873
10	0.060	0.497239
14	0.065	0.497279
9	0.070	0.497182
15	0.075	0.496999
8	0.080	0.497642
16	0.085	0.496556
3	0.095	0.495844
4	0.096	0.496167
5	0.097	0.496018
6	0.098	0.496273
7	0.099	0.496013
0	0.110	0.495383
1	0.120	0.495062
2	0.150	0.492045